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1 .  ESTIMATION  METHODS 


The  position  of  a  stationary  transmitter  or  radiating  emitter  can  be 
estimated  from  passive  measurements  of  the  arrival  times,  directions  of  arriv¬ 
al,  or  Doppler  shifts  of  electromagnetic  waves  received  at  various 
stations.1'2  In  the  first  three  sections,  the  basic  methods  of  estimation 
applicable  to  transmitter  location  and  the  accuracy  of  estimators  are  exam¬ 
ined.  Passive  location  systems  using  arrival-time  and  bearing  measurements 
are  considered  in  detail  in  sections  4  and  5,  respectively.  The  use  of 
Doppler  information  is  briefly  summarized  in  section  6. 

The  components  of  an  n-dimensional  vector  x  that  is  to  be  estimated  are 
the  position  coordinates  in  two  or  three  dimensions  and  possibly  other  param¬ 
eters  such  as  the  time  of  emission  of  the  radiation.  A  set  of  N  measurements 
r^,  i  1 ,  2,  .  .  .,  N,  is  collected  at  various  positions.  In  the  absence  of 
random  measurement  errors,  r^  is  equal  to  a  known  function  f^(x).  In  the 
presence  of  additive  errors, 

ri  ”  *  ni  ,  i“1,  2,  ...,  N  .  (1) 

These  N  equations  can  be  written  as  a  single  equation  for  N-dimensional  column 
vectors: 


r  -  f{«)  +  n  .  (2) 

The  measurement  error  n  is  assumed  to  be  a  multivariate  random  vector  with  an 
N  x  n  positive-definite  covariance  matrix 

n  -  E[(n  -  E[n])(n  -  E[n])T]  ,  (3) 

where  e[  ]  denotes  the  expected  value  and  the  superscript  T  denotes  the  trans¬ 
pose. 

If  x  is  regarded  as  an  unknown  but  nonrandom  vector  and  n  is  assumed  to 
have  a  zero  mean  and  a  Gaussian  distribution,  then  the  conditional  density 
function  of  r  given  x  is 

pdrlx)  -  -  "-n7^|--|-i7^  exp  {-  i  [r  -  f(x)]TH~1[r  -  f  («)  ]}  ,  (4) 

where  |h|  denotes  the  determinant  of  ■,  and  the  superscript  -1  denotes  the 
inverse.  Because  M  is  symmetric  and  positive  definite,  its  inverse  exists. 
The  maximum-likelihood  estimator  is  that  value  of  x  which  maximizes  equation 
(4).  Thus,  the  maximum- likelihood  estimator  minimizes  the  quadratic  form 


1L.  H.  Wegner,  On  the  Accuracy  Analysis  of  Airborne  Techniques  for  Passive¬ 
ly  Locating  Electromagnetic  Emitters,  Rand  Corp .  R-722-PR,  Nat,  Tech.  Inf. 
Serv.  AD  729  767  (1971). 

2ff.  B.  Lee,  A  Nor  el  Procedure  for  Assessing  the  Accuracy  of  Hyperbolic 
Mu ltllateratlon  Systems,  IEEE  Trans.  Aerosp.  Electron.  Syst.  AES- 11  (January 
1975),  2. 
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(5) 


Q(x)  =  [r  -  f(x)]TN_1[r  -  f (x)  j 


Minimization  of  Q(x)  is  a  reasonable  criterion  for  determination  of  an  esti¬ 
mator  even  when  the  additive  error  cannot  be  assumed  to  be  Gaussian.  In  this 
case,  the  resulting  estimator  is  called  the  least-squares  estimator  and  is 
regarded  as  a  matrix  of  weighting  coefficients. 

In  general,  f(x)  is  a  nonlinear  vector  function.  To  determine  a  reason¬ 
ably  simple  estimator,  f(x)  can  be  linearized  by  expanding  it  in  a  Taylor 
series  about  a  reference  point  specified  by  the  vector  Xq  and  retaining  the 
first  two  terms;  that  is,  we  use 


£(x)  =  f(xo)  +  G(x  -  Xq) 


(6) 


where  x  and  Xq  are  n  x  i  column  vectors  and  G  is  the  N  x  n  matrix  of  deriva¬ 
tives  evaluated  at  Xgi 


3fN 

X-Xq 


ill 

**n 


3x_ 


X4r 


(7) 


Each  row  of  this  matrix  is  the  gradient  vector  of  one  of  the  components  of 
f(x).  The  vector  Xg  could  be  an  estimate  of  x  determined  from  a  previous 
iteration  of  the  estimation  procedure  or  based  upon  a  priori  information.  It 
is  assumed  in  the  subsequent  analysis  that  Xg  is  sufficiently  close  to  x  that 
equation  (6)  is  an  accurate  approximation. 


Combining  equations  (5)  and  (6)  gives 

Q(x)  =  (r,  -  Gx)TN"1(r1  -  Gx)  , 

where 

r1  -  r  ~  f(*o)  +  <**0  * 


(8) 


(9) 


To  determine  the  necessary  condition  for  the  estimator  x  that  ^minimizes  Q(x), 
we  calculate  the  gradient  of  Q(x),  defined  by 


T  3Q 

7xQ{x)  = 


3Q  3Q  "It 

177  •  •  •  JTn\  ' 


(10) 


and  then  solve  for  the  x  such  that  7  Q(x)  ■  0.  From  its  definition. 


symmetric  matrix;  that  is,  NT  =  N.  Since  (h  * )  x  t  _ 

(H'lf  =  N_1,  which  implies  that  H” 1  is  a  symmetric  matrix.  Therefore, 


?xQ<*>!x=£ 


0 


(11) 


»  2GTIT 1 Gx  -  2GTH“lr1  - 

We  assume  that  the  matrix  GTH_1G  is  nonsingular.  Thus,  the  solution  of  equa¬ 
tion  (11)  is 

;  =  (G'VlGpG'Vlr, 

=  Xq  +  (gtn“1g)“1GTN-1 [r  -  f(xQ)]  .  (12) 

Using  equation  (12),  direct  calculation  shows  that  equation  (8)  can  be  written 
in  the  form 

Q(x)  =  [x  -  ijTGTN-lG[x  -  x\  -  r>-lG(GTM-lG)'lri  +  ,  (13) 

where  only  the  first  term  depends  upon  x.  Since  H  is  symmetric  and  positive 
definite,  it  has  positive  eigenvalues.  If  He  -  Xe,  then  N_1e  *  X_1e.  Thus, 
if  e  is  an  eigenvector  of  N  with  eigenvalue  X,  then  e  is  also  an  eigenvector 
of  with  eigenvalue  1/X.  Since  it  is  symmetric  and  its  eigenvalues  are 

positive,  H_1  is  positive  definite.  Therefore,  x  =  x  minimizes  Q(x).  The 
estimator  of  equation  (12)  is  called  the  linearized  least-squares  estimator. 


Substituting  equation  (2)  into  equation  (12)  and  rearranging  terms,  the 
expression  for  x  can  be  written  in  the  form 

x  =  x  +  (gtH_1g)"1Gtw"1 [f (x)  -  f (xq)  -  G(x  -  Xq)  +  n]  ,  (14) 

which  shows  how  the  estimator  error  is  affected  by  the  linearization  error  and 
the  noise.  The  bias  of  the  estimator  x  is  defined  as  b  ■  e[x]  -  x.  Using 
equation  (14),  we  obtain 

b  -  (GTB"1G)_1GTH_1{f(*)  ~  f(*o)  "  G(x  "  xo)  +  Elnlf  ’  (15) 

If  equation  (6)  is  exact  and  E[nl  »  0,  then  the  least-squares  estimator  is 
unbiased.  If  systematic  errors  occur  in  the  measurements,  then  Eln]  *  0.  To 
minimize  the  estimator  bias  due  to  systematic  errors,  the  magnitude  of  each 
E[njJ  should  be  minimized  through  system  calibrations.  If  some  of  the  E[n^] 
are  known  functions  of  various  parameters  and  N  is  sufficiently  large,  these 
parameters  can  be  made  components  of  the  vector  x  and  estimated  along  with  the 
other  components  of  x.  The  bias  due  to  the  nonlinearity  of  f(x)  can  be  esti¬ 
mated  by  expanding  f(x)  in  a  Taylor  series  about  Xq  and  retaining  second-order 
terms. 


can  compute  both  estimate  and  covariance  simultaneously.  If  n  is  zero-mean 
Gaussian,  the  maximum- likelihood  or  least-squares  estimator  for  the  linearized 
model  is  the  same  as  the  minimum-variance  unbiased  estimator.3 4 

The  measurement  error  vector  n  is  assumed  to  encompass  all  the  contribu¬ 
tions  to  error,  including  uncertainties  in  the  system  or  physical  parameters, 
such  as  the  station  coordinates  or  the  speed  of  propagation.  If  q  is  a  vector 
of  the  parameters,  then  the  measurement  vector  r  can  often  be  expressed  as 

r  ■  f1 (x,q)  +  n1  ,  (17) 

where  f  1  (  )  is  a  vector  function  and  n.  is  the  random  error  due  to  causes 

unrelated  to  uncertainties  in  q.  Let  qQ  denote  the  assumed  value  of  q.  If 
is  sufficiently  close  to  q,  then  a  Taylor  series  expansion  yields 

f^x.q)  =  f i (x* q<) )  +  Gl(q  "  **0^  '  (18) 

where  G1  is  the  matrix  of  derivatives  with  respect  to  q  evaluated  at  qQ. 
Equation  (2)  results  from  making  the  identifications 

f(x)  =  f.j (x»qg)  ,  n  =  Gt(q  -  qQ)  +  n,  .  (19) 

If  q  is  nonrandom,  then  the  parameter  uncertainties  ultimately  contribute  to 
the  bias  of  the  least-squares  estimator.  If  q  is  random,  then  the  variance 
and  possibly  the  bias  are  affected. 

Any  a  priori  information  can  be  incorporated  into  the  estimation  procedure 
in  several  ways.  It  can  be  used  to  select  an  accurate  reference  point  Xq  for 
the  first  iteration  of  the  least-squares  estimator.  If  the  transmitter  is 
known  to  be  located  within  a  region,  but  the  estimated  position  is  outside 
this  region,  a  logical  procedure  is  to  change  the  estimate  to  the  point  in  the 
region  that  is  closest  to  the  original  estimate.  If  an  a  priori  dittribution 
function  for  the  transmitter  position  can  be  specified,  a  Bayesian  estimator 
can  be  determined.  However,  the  Bayesian  estimator  is  usually  too  complex  a 
mathematical  function  to  yield  a  feasible  computational  algorithm  unless 
simplifying  assumptions  are  made  about  the  a  priori  distribution.1* 

The  location  estimate  can  be  continually  refined  if  a  sequence  of  measure¬ 
ments  is  taken.  If  successive  measurements  are  uncorrelated,  a  new  least- 
squares  estimate  can  be  determined  by  combining  new  measurements  with  the  old 
estimate.3  Since  measurements  do  not  have  to  be  stored  after  processing,  a 
significant  computational  savings  is  sometimes  possible. 


3A.  P.  Sage  and  J.  L.  Melsa,  Estimation  Theory  with  Applications  to  Commu¬ 
nications  and  Control,  McGraw-Hill  (1971). 

4P.  J .  Butterly,  Position  Finding  with  Empirical  Prior  Knowledge,  IEEE 
Trans.  Aerosp.  Electron.  Syst.  AES-8  (March  1972),  142. 


2 .  ESTIMATOR  ACCURACY 


If  r  is  a  Gaussian  random  vector,  then  equation  (12)  indicates  that  x  is  a 
Gaussian  random  vector.  Its  probability  density  function  is 

f£<5>  =  [(  2tt  )  n/2  |  p|  l/2j- 1  exp  [-1  (c  -  m)Tp-1(c  -  ■)]  ,  (20) 

where  ■  =  E[x]  is  the  mean  vector,  and 

P  =  E[(x  -  ■) (x  -  -)T]  (21 ) 

is  the  covariance  matrix  given  by  equation  (16).  By  definition,  P  is  symmet¬ 
ric  and  positive  semidef inite .  Thus,  it  has  nonnegative  eigenvalues.  Equa¬ 
tion  (16)  indicates  that  P  *  exists.  Therefore,  P  does  not  have  zero  as  an 
eigenvalue.  Thus,  P  is  positive  definite. 

The  loci  of  constant  density  function  values  are  described  by  equations  of 
the  form 


(e  -  «)Tp-1(e  -  ■)  =  <  ,  (22) 

where  k  is  a  constant  that  determines  the  size  of  the  n-dimensional  region 
enclosed  by  the  surface.  In  two  dimensions,  the  surface  is  an  ellipse;  in 
three  dimensions,  it  is  an  ellipsoid;  in  the  general  case  of  n  dimensions,  it 
may  be  considered  a  hyperellipsoid.  Unless  P  is  a  diagonal  matrix,  the  prin¬ 
cipal  axes  of  the  hyperellipsoids  are  not  aligned  with  the  coordinate  axes. 

A 

The  probability  that  x  lies  inside  the  hyperellipsoid  of  equation  (22)  is 

Pe(tc)  =//.../  f“(£)  d£i  d^2  •  •  •  dCn  ,  (23) 

R 

where  the  region  of  integration  is 

R  =  U  :  <5  *  ■>TP"i<5  -■><<}  •  (24) 

To  reduce  equation  (23)  to  a  single  integral,  we  perform  a  succession  of 
coordinate  transformations.  First,  we  translate  the  coordinate  system  so  that 
its  origin  coincides  with  m  by  making  the  change  of  variables  y  =  £  - 
Since  the  Jacobian  is  unity,  we  obtain 

Pe(K>  =  a  //  .  .  .  /  exp(-  4j  YTP-1y)  dy,  <*Y2  .  .  .  dYn  ,  (25) 

R1 

where 


{y  :  YTP_1Y  <  <}  r 

(26) 

[(2x)n/2|p|l/2]-l  . 

(27) 

To  simplify  equation  (25),  we  rotate  the  coordinate  axes  so  that  they  are 
aligned  with  the  principal  axes  of  the  hyperellipsoid.  Because  P  is  a  symmet¬ 
ric  positive-definite  matrix,  so  is  P_1.  Therefore,  an  orthogonal  matrix  A 
(with  eigenvectors  as  columns)  exists  that  diagonalizes  P_1.  Thus,  AT  *  A-1 


r- 


r 


! 


atp-1a  = 


x^1 


x"1 


-  [X"1]  . 


(28) 


where  X1#  X2,  .  .  Xn  are  the  eigenvalues  of  P.  A  rotation  of  axes  results 

in  new  variables  defined  by 


T 

C  =  AY 


(29) 


Since  ATA  =  I  and  the  determinant  of  the  product  of  matrices  is  equal  to  the 
product  of  the  determinants  of  the  matrices,  the  determinant  of  AT,  which  is 
the  Jacobian  of  the  transformation,  is  unity.  Substituting  equations  (28)  and 
(29)  into  equations  (25)  and  (26)  yields 


Pe(*c)  =  a  ff  ...  f  exp(-  CT[X-1]c)  d?,  d?2  .  .  .  d^n 
R2 

,  /  exp^-  \  d^2  *  *  *  d^n  ' 


(30) 


=  a 


II . 


where 


j.sf 


S  K 


(31) 


and  the  ^  are  the  components  of  Region  R2  is  the  interior  of  a  hyper¬ 

ellipsoid  with  principal  axes  of  lengths  2/ kX^ ,  i  =  1,  2,  .  .  .,  n.  By  intro¬ 
ducing  new  variables 


where  T(  )  is  the  gamma  function.  Therefore,  the  differer 
p  and  p  +  dp  is 


volume  between 


and  equation  (33)  can  be  reduced  to 


where  the  error  function  is  defined  by 


Equation  (40)  is  obtained  by  integration  by  parts 


To  verify  equation  (35),  we  define  the  volume  of  a  hypersphere  of  radius  p 


A  change  of  coordinates  shows  that 


where  vnO)  is  the  volume  of  a  unit  hypersphere.  Straightforward  calculations 
give 


v^l)  -  2 


v2(i)  -  * 


(44) 


where  the  "volumes"  are  a  length  and  an  area,  respectively.  We  define  the 
sets 


1(X1'X2)  :  X1  +  X1  ^  ’1  ' 

n 

(x3,  .  .  xn)  :  y  x|  i  1 
1=3 


-  *2 


(45) 

(46) 


For  n  >  3,  Fubini's  theorem  for  interchanging  the  order  of  integrations  and  a 
change  of  coordinates  in  equation  (42)  give 


Vn( 1 )  -  /  /  dx1  dx2  //  •  •  •  /  dx3  ...  dxn 

c(x1 ,x2) 

Vn-2  (VTTTrrri)  dx,  dx2  . 
Equation  (43)  and  further  coordinate  changes  yield 

Vn(1)  "  Vn-2(1)  f  A1  *  X1  -  x£)(n~2)/2  <*xi 

=  vn_2(i)  rr  (l  -  r2)(n-2)/2  r  dr  d0 

r 

Jo 


»  i  (n-2)/2 
xV  (1)  |  x  dx 

n-z 


2irV  .(1) 
n-z 


By  induction,  this  recursion  relation  and  equations  (44)  imply  that 


2L2ll 


m-1 


v2m(1)  *  2  •  4  .  .  .  (2m)  •  v2m-1(1)  *  1  •  3  ...  (2m  -  1)  » 

m  *  1 ,  2,  .  .  . 


(47) 


(48) 


(49) 


We  can  express  Vn(1)  in  terms  of  a  compact  formula  by  using  the  properties  of 
the  gamma  function:  T(t  +  1)  ■  tT(t)>  T(1)  »  1;  1*0/2)  *  /»•  We  obtain 


,"/2 


Vn(1)  "  r(*  ♦  i) 


,  n  ■  1,  2, 


(50) 
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Combining  equations  (43)  and  (50)  yields  equation  (35). 

If  Pg  is  specified,  say  Pg  =  1/2,  then  equation  (37),  (38),  (39),  or  (40) 
can  be  solved  numerically  to  determine  the  corresponding  value  of  x,  which  in 
turn  defines  a  hyperellipsoid  by  equation  (22).  The  concentration  ellipsoid 
corresponding  to  probability  Pg  is  defined  to  be  the  particular  hyperellipsoid 
for  which  Pg  is  the  probability  that  x  lies  inside  it.  Thus,  the  concentra¬ 
tion  ellipsoid  is  a  multidimensional  measure  of  accuracy  for  an  unbiased 
estimator. 

A  scalar  measure  of  estimator  accuracy  is  the  root-mean-square  error,  er, 
which  is  defined  by 

er  =  E[.|1  (*i  *  xi)2]  *  (51) 

Expanding  equation  (51)  and  using  equation  (21),  we  obtain 

n 

e l  =  tr (P)  +  I  b?  ,  (52) 

r  i  =  1  1 

where  tr(P)  denotes  the  trace  of  P  and  b^  =  e[x^]  -  x^  denotes  a  component  of 
the  bias  vector  b. 


3.  TWO-DIMENSIONAL  ESTIMATORS 


For  the  estimator  of  a  two-dimensional  vector,  such  as  position  coordi¬ 
nates  on  the  surface  of  the  earth,  the  bivariate  covariance  matrix  can  be  ex¬ 
pressed  as 


P  = 


L°12 


(53) 


A  straightforward- calculation  yields  the  eigenvalues: 

*1  =  [of  +  o|  +  V(°1  "  °2 +  4o12^J  '  (54> 

*2  "  2  [°1  +  °2  ~>l(a  1  ~  °i)2  +  4o1 2]  '  (55) 

where  the  positive  square  root  is  used.  By  definition,  A1  >  X2* 

Suppose  that  new  coordinates  are  defined  by  rotating  the  axes  of  the  old 

coordinate  system  counterclockwise  through  an  angle  0  as  shown  in  figure  1 .  A 
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vector  represented  by  y  in  the  old  coordinates  is  represented  in  the  new 
coordinates  by  C  =  ATy,  where  A  is  the  orthogonal  matrix 


Figure  1 .  Concentration  ellipse  and 
coordinate  axes. 


From  equations  (53)  and  (56),  direct  but  lengthy  calculation  shows  that  ATP_1A 
is  a  diagonal  matrix  and  the  columns  of  A  are  eigenvectors  if 


-  i  <  6  i 


(57) 


If  =  o|  and  =  0,  we  take  8  =  0,  Since  the  determinant  of  a  matrix  is 
equal  to  the  product  of  the  eigenvalues,  =  a*o|  -  o^.  Using  this 

result,  the  diagonal  matrix  can  be  written  in  the  form 


[V  0  "I 

[x  1 ]  =  ,  f  (58) 

o  a-i_ 


Since  P"1  exists  according  to  equation  (16),  neither  eigenvalue  can  equal  zero 
and  [X_1]  is  well  defined. 


A  concentration  ellipse  defined  by  yAP“1Y  =  k  in  the  old  coordinates  is 
described  by  (l^/Xj)2  +  ( ^2 )  *  K  or  (^1/^2)2  +  (Cj/^^2  =  k  in  the  new 
coordinates,  a  fact  which  indicates  that  the  new  axes  coincide  with  the  prin¬ 
cipal  axes  of  the  ellipse.  Thus,  equation  (57)  represents  the  angular  offset 
of  one  of  the  principal  axes  of  the  ellipse  relative  to  the  old  coordinate 
axes.  Figure  1  depicts  a  concentration  ellipse  and  the  appropriate  angle  of 
axis  rotation.  Since  X1  >  X2,  the  major  and  minor  axes  have  lengths  2/kK 1  and 
2/ kX2»  respectively.  If  the  ellipse  encloses  a  region  that  includes  a 
Gaussian  random  vector  with  probability  Pfi,  then  equation  (39)  implies  that 

*c  =  -2Ln(l  -  Pft)  .  (60) 

Suppose  that  a  two-dimensional  Gaussian  random  vector  describes  the  esti¬ 
mated  location  of  a  transmitter.  A  crude  but  simple  measure  of  accuracy  is 
the  circular  error  probable  (CEP)  (often  erroneously  called  the  circular  error 
probability,  despite  the  fact  that  it  is  not  a  probability).  The  CEP  is 
defined  as  the  radius  of  the  circle  that  has  its  center  at  the  mean  and  con¬ 
tains  half  the  realizations  of  the  random  vector.  The  CEP  is  a  measure  of  the 
uncertainty  in  the  location  estimator,  x,  relative  to  its  mean,  e[x].  If  the 
location  estimator  is  unbiased,  the  CEP  is  a  measure  of  the  estimator  uncer¬ 
tainty  relative  to  the  true  transmitter  position.  If  the  magnitude  of  the 
bias  vector  is  bounded  by  B,  then  with  a  probability  of  one  half,  a  particular 
estimate  is  within  a  distance  of  B  +  CEP  from  the  true  position.  The  geomet¬ 
rical  relations  are  depicted  in  figure  2. 


PARTICULAR 

ESTIMATE 


Figure  2.  Geometry  of  transmitter  position,  mean 
location  estimate,  CEP,  estimator  bias  vector,  and 
a  particular  location  estimate. 


From  the  definition,  it  follows  that  we  can  determine  the  CEP  by  solving 
the  equation 

4  -  /  /  fJU>  <3?i  <352  »  (61) 

R 

where 
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R  -  {t  *  U  -  »|  <  cep}  . 


In  a  manner  analogous  to  the  derivati  'n  of  equation  (30),  we  successively 
translate  and  rotate  coordinates  to  obtain 


where 


*  ’  1^/Ri/“P('  *  ji  $ )  d{’  dc>  ' 


R1  -  {(?1»C2)  :  (e?  +  c|)l/2  <  CEP}  , 


and  the  are  given  by  equations  (54)  and  (55).  Changing  to  polar  coordi¬ 
nates  by  substituting  51  *  r  cos  0  and  «;2  «  r  sin  0,  we  get 


■  i:t  r  ™[-  £  Hr*  *  ^)] 


dr  d0  . 


To  simplify  equation  (65),  we  do  some  preliminary  manipulations.  ‘Ole 
modified  Bessel  function  of  the  first  kind  and  zero  order  can  be  expressed  as 


Ig (x)  “'2 H  Jq  exP(x  cos  e)  d0 


Because  of  the  periodicity  of  the  integrand,  we  also  have 

.  f 2x(n+1) 

I0(x)  »  /  exp(x  cos  0)  d0 

2x n 


for  any  integer  n.  Adding  m  equations  of  this  form  with  successive  values  of 
n,  we  obtain 


.  f  2irm 

mlg(x)  *  J  exp(x  cos  0)  d0  , 


m  =  1 ,  2, 


Changing  coordinates  with  0  =  m$  gives 


Ig ( x )  *  ^  exp(x  cos  m$)  d$  ,  m  »  1,  2,  . 


Trigonometric  identities  yield 


ssarf- +  -  2^7 +  2^  +  (2I7  -  2^) cos  26 


(70) 


Substituting  equation  (70)  into  equation  (65)  and  using  equation  (69),  we 
obtain 


^  -  Jo  r  exp  -{A *  ’<>[(^1 '  ^r)r2J 


dr  . 


(71) 


A  final  change  of  coordinates  yields 


p 

-4  (1  ♦  Y2)  -  / 

4y2  Jo 


+Y2) 


exp(-x)  I 


9  ^2 

dx,  Y2  =  i 
A1 


(72) 


The  form  of  this  relation  implies  that  the  CEP  has  the  form  CEP  *  /X^  f(y)  for 
some  function  f(  ).  If  a12  *  0  an<*  °i  =  °2  =  a’  then  ^1  =  ^2  *  °2  an<*  e1ua~ 
tion  (72)  can  be  solved  to  show  that  CEP  *■  1.177o.  In  the  general  case  where 
X1  *  X2,  numerical  integration  is  necessary  to  solve  for  the  CEP.  A  simple 
approximation  is 


CEP  -  0.563-^7^  +  0.614-^XJ 


(73) 


which  is  accurate  to  within  1  percent  for  y  “  0.3  or  larger,  underestimates 
the  CEP  by  less  than  10  percent  for  0.1  <  Y  <  0.3,  and  underestimates  by  less 
than  20  percent  elsewhere.  Although  approximations  that  are  more  accurate  for 
small  Y  are  easily  produced,  they  are  usually  irrelevant  because  the  eccentri¬ 
city  of  the  concentration  ellipse  for  small  y  may  be  too  pronounced  for  the 
CEP  to  be  an  adequate  performance  measure.  An  approximation  that  is  accurate 
to  within  approximately  10  percent  for  all  values  of  y  is 

CEP  «  0.75-y^  +  *2  =  °*75V°i  +  ,  (74) 

where  the  last  relation  follows  from  the  fact  that  the  trace  of  a  matrix  is 
equal  to  the  sum  of  its  eigenvalues.  Above  y  »  0.4,  this  approximation  under¬ 
estimates  the  CEP;  below  y  "  0.4,  it  overestimates  the  CEP.  For  an  unbiased 
estimator,  equation  (52)  implies  that  CEP  «  0.75  er . 

4.  HYPERBOLIC  LOCATION  SYSTEMS 

Hyperbolic  location  systems,  which  are  often  called  TDOA  (time  difference 
of  arrival)  systems,  locate  a  transmitter  by  processing  signal  arrival-time 
measurements  at  three  or  more  stations.  Measurements  at  two  stations  are 
combined  to  produce  a  relative  arrival  time  that,  in  the  absence  of  noise  and 
interference,  restricts  the  possible  transmitter  location  to  a  hyperboloid 
with  the  two  stations  as  foci.  Transmitter  location  is  estimated  from  the 
intersections  of  three  or  more  independently  generated  hyperboloids  determined 
from  at  least  four  stations.  If  the  transmitter  and  the  stations  lie  in  the 


same  plane,  location  is  estimated  from  the  intersections  of  two  or  more  hyper¬ 
bolas  determined  from  at  least  three  stations,  as  illustrated  in  figure  3. 
Two  hyperbolas  may  have  more  than  one  point  of  intersection.  The  resulting 
location  ambiguity  may  be  resolved  by  a  priori  information,  bearing  measure¬ 
ments  at  the  stations,  or  the  use  of  additional  hyperbolas. 


HYPERBOLA  FROM  HYPERBOLA  FROM 

STATION  1  AND  2  STATION  2  AND  3 


Figure  3.  Intersecting  hyperbolas  from  three 
stations . 

Suppose  that  the  arrival  times  t1 ,  t2»  .  .  .,  t^  of  a  signal  transmitted 
at  time  tQ  are  measured  at  N  stations  having  positions  specified  by  the  column 
vectors  ,  s2,  .  .  .,  Sjj.  The  geometrical  configuration  is  illustrated  in 

figure  4.  The  stations  may  be  at  fixed  positions  or  represent  various  points 
along  the  trajectory  of  an  aircraft.  If  the  signal  velocity  is  c  and  is 
the  propagation  path  length  between  the  transmitter  and  station  i,  then 

Di 

fci  "  fc0  +  c  +  ei  *  i  -  1,  2,  .  .  .,  N  .  (75) 


TRANSMITTER 


Figure  4.  Geometry  of  transmitter 
and  N  stations. 
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The  arrival-time  measurement  error  accounts  for  propagation  anomalies, 
receiver  noise,  and  errors  in  the  assumed  station  positions.  In  matrix  form, 
equation  (75)  becomes 

t  =  tQ1  +  £  D  +  e  ,  (76) 

where  t,  D,  and  e  are  N-dimensional  column  vectors  with  components  t^,  D^, 
and  e. ,  i  =  1 ,  2,  .  .  .,  N,  respectively,  and  1  is  a  column  vector  of  ones. 
Suppose  that  we  seek  to  estimate  both  tg  and  the  column  vector  R,  with  compo¬ 
nents  x,  y,  and  z,  that  specifies  the  transmitter  position.  Equation  (76)  has 
the  form  of  equation  (2)  with  r  =  t,  f(x)  =  tQ1  +  D/c,  n  =  e,  and  x  * 
[tp  x  y  z]^.  for  line-of-sight  propagation  from  the  transmitter  to  the  sta¬ 
tions,  Di  =  |R  "  sil*  where  I  I  represents  the  Euclidean  norm.  Let  the 
column  vector  R^,  with  components  xQ,  yg,  and  zQ,  specify  a  reference  point 
near  the  transmitter  position.  Let  Dg^  =  |Rg  -  s^J  denote  the  distance  from 
station  i  to  the  reference  point.  Using  equation  (7)  with  Xg  =  [o  Xg  yg  Zg]T 
after  expressing  each  JR  -  s^|  in  terms  of  its  components,  we  obtain 

G  =  [ 1  P/c]  ,  (77) 


where 


(*o  "  »i)T/°oi 

(*0  "  *n)T/DON 


(78) 


Each  row  of  P  is  the  unit  vector  pointing  from  one  of  the  stations  to  the 
reference  point.  Equation  (12)  with  the  above  relations  and  substitutions 
gives  the  least-squares  or  maximum-likelihood  estimator;  equation  (16)  pro¬ 
vides  the  covariance  matrix  of  the  estimator. 

In  hyperbolic  systems,  no  attempt  is  made  to  estimate  tg.  We  eliminate  it 
from  consideration  by  measuring  the  relative  arrival  times: 


ci+1 


(°i  -  Pl+J 


+  n. 


2,  .  .  .,  N  -  1 


(79) 


where  n^  is  the  measurement  error.  Measuring  time  differences  is  not  the  only 
way  to  eliminate  tg,  but  it  is  the  simplest.  If  the  relative  arrival  times 
are  determined  by  subtracting  measured  arrival  times,  then 

n^  ■  —  e^^  ,  i  ■  1 ,  2,  .  .  .,  N  —  1  .  (80) 

The  n^  have  zero  means  if  successive  have  equal  means,  even  if  the  latter 
means  are  nonzero.  A  nonzero  E[n^J  may  result  from  uncalibrated  different 
time  delays  or  unsynchronized  clocks  in  two  receivers.  If  the  relative  arriv¬ 
al  times  are  determined  by  cross  correlation,  then  equation  (80)  is  not  neces¬ 
sarily  valid. 
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If  the  transmitter  produces  a  sequence  of  pulses,  the  corresponding  re¬ 
ceived  pulses  at  stations  i  and  i  +  1  must  be  correctly  associated  in  measur¬ 
ing  the  time  difference  t^  -  t^+1 .  A  potential  ambiguity  arises  when  the  time 
difference  exceeds  the  time  between  successive  pulse  transmissions.  This 
ambiguity  may  be  resolved  by  using  bearing  measurements  or  a  priori  informa¬ 
tion  to  eliminate  associations  that  lead  to  impossible  location  estimates. 

In  matrix  form,  equation  (79)  becomes 

Ht  =  ~  ID  +  n  ,  (81 ) 

where  we  use  the  (N  -  1 )  x  N  matrix 


H  = 


1  -1  0 
0  1  -1 


_0  0  0 
If  equation  (80)  is  valid,  then 


(82) 


n  =  He  .  (83) 

Since  we  seek  to  estimate  the  position  vector  R,  equation  (81)  has  the 
form  of  equation  (2)  with  r  =*  Ht,  f(x)  *  HD/c,  and  x  =  R.  A  direct  calcula¬ 
tion  of  G  yields 

G  =  £  HF  ,  (84) 


where  F  is  defined  by  equation  (78).  Let  N£  denote  the  covariance  matrix  of 
the  arrival-time  errors.  If  equation  (83)  holds,  then  the  covariance  matrix 
of  the  measurement  errors,  defined  by  equation  (3),  is  related  to  H£  by 

H  =  HNeHT  .  (85) 

Using  equation  (84),  equation  (12)  implies  that  the  least-squares  estimator  is 

R  =>  Rq  +  c(FTHTN"1HF)*1FTHTH-1(Ht  -  £  HDQ)  ,  (86) 

where  DQ  has  components  D0i'  =  1  *  2,  •  •  •»  N*  The  estimator  is  unbiased  if 
n  is  a  zero-mean  random  variable  and  the  linearization  error  is  negligible. 
The  covariance  matrix  of  R,  given  by  equation  (16),  is 

P  =  c2(FTHTII~1HF)"1  .  (87) 

Equation  (86)  is  valid  for  line-of-sight  propagation.  If  the  signal 
propagation  to  the  stations  involves  atmospheric  reflections,  the  equations 
for  the  change  and  thus  the  estimator  changes. 
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In  general,  the  least-squares  estimator  requires  knowledge  of  the  statis¬ 
tics  of  the  measurement  errors.  However,  if  equation  (85)  applies,  if  the 
covariances  of  the  are  zero,  and  if  the  variances  of  the  have  the  common 
value  a2,  then  cancellation  in  equation  (86)  leaves  an  estimator  that  is 
independent  of  o2.  Equality  of  the  variances  is  a  reasonable  assumption  for 
stations  with  identical  receivers  that  are  much  closer  to  each  other  than  to 
the  transmitter. 

Let  denote  the  variance  of  the  measured  arrival  time  t^  at  station 
i.  When  equation  (83)  is  valid,  the  mean-square  ranging  error  is  defined 
as  c2o2,  where 

°i  =  i  1,  °li  <88) 

is  the  average  variance  of  the  arrival  times.  The  geometric  dilution  of 
precision  (GDOP)  is  defined  as  the  ratio  of  the  root-mean-square  position 
error,  ef,  to  the  root-mean-square  ranging  error.  It  follows  from  equation 
(52)  that  the  GDOP  associated  with  an  unbiased  estimator  and  a  hyperbolic 
system  is 


GDOP  = 


'trace  I 
co„ 


The  GDOP  indicates  how  much  the  fundamental  ranging  error,  intuitively  meas¬ 
ured  by  cog,  is  magnified  by  the  geometric  relation  among  the  transmitter 
position  and  the  stations.  If  the  geometry  is  such  that  the  arrival-time 
variances  are  nearly  equal,  then  the  GDOP  is  only  weakly  dependent  on  them. 
For  the  two-dimensional  location  problem,  equations  (74)  and  (89)  yield 

CEP  «  (0.75cog)GDOP  .  (90) 

Since  the  arrival- time  variance  o|  is  due  primarily  to  the  thermal  and 
environmental  noise,  it  is  often  reasonable  to  model  as  the  sum  of  a  con¬ 
stant  bias  plus  zero-mean  white  Gaussian  noise.  The  Cramer-Rao  bound  for  an 
arrival-time  estimate  in  the  presence  of  white  Gaussian  noise  gives^ 


2  (f  B?)M 


where  E  is  the  energy  in  the  received  signal,  NQ/2  is  the  two-sided  noise 
power  spectral  density,  and  82  is  a  function  of  the  bandwidth  of  the  signal. 
If  S(o))  denotes  the  Fourier  transform  of  the  signal,  then 

,  /*„,  U)2  I  S  (  U)  )  I  2  dcu 

e2  -  —z  7 — r; -  •  (92) 

/_»  l  S  ( a) )  |  2  dto 


If  the  received  signal  consists  of  pulses,  then  E  is  the  sum  of  the  energies 
of  the  individual  pulses.  An  approximate  model  for  many  radar  signals  is  a 
series  of  pulses,  each  of  which  results  from  passing  a  truncated  sinusoid  with 


A.  D.  Whalen,  Detection  of  Signals  in  Noise,  Academic  Press  (1971) 
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an  ideal  rectangular  envelope  of  duration  Tp  through  an  ideal  rectangular 
bandpass  filter  of  bandwidth  B  centered  at  the  sinusoidal  frequency.  For  each 
pulse  and  for  the  entire  radar  signal,  equation  (92)  yields 

02  „  2B/Tp  ,  BTp  »  1  .  (93) 

In  contrast,  for  a  signal  with  a  uniform  Fourier  transform  over  a  bandwidth  B, 
equation  (92)  gives 

02  =  7r2B2/3  .  (94) 

This  model  might  approximate  a  communications  signal. 

Let  T  denote  the  total  signal  duration,  Rg  =  E/T  denote  the  average  signal 
power  at  the  receiver,  and  D  denote  the  distance  between  the  transmitter  and 
the  receiver.  Over  a  large  range  of  values  of  D,  it  is  often  possible  to 
approximate  Rg  by6 

K  exp(-aD) 

R  =  — -  ,  (95) 

3  D" 

where  a,  n,  and  KE  are  independent  of  D,  but  may  be  functions  of  other  param¬ 
eters  such  as  the  transmitter  power,  antenna  gains,  antenna  heights,  and  the 
signal  frequency.  For  optical  and  millimeter-wave  frequencies,  accurate 
modeling  requires  a  >  0,  but  we  may  usually  set  a  =  0  at  other  frequencies. 
Inequality  (91)  and  equation  (95)  relate  o£  to  D. 

As  an  important  special  case,  we  consider  a  transmitter  and  three  stations 
in  the  same  plane  so  that  only  two  position  coordinates  are  to  be  estimated. 
The  planar  model  is  reasonable  if  a  transmitter  and  stations  are  near  the 
surface  of  the  earth  and  close  enough  that  the  curvature  of  the  earth's  sur¬ 
face  can  be  neglected.  One  of  the  stations  is  designated  the  master  station, 
and  the  other  two  are  called  slave  stations.  Arrival-time  measurements  at  the 
slave  stations  are  sent  to  the  master  station,  where  the  time  differences  and 
then  the  position  estimate  are  computed. 


We  assume  that  the  are  uncorrelated  random  variables  so  that 


6D.  J,  Torrieri,  Principles  of  Military  Communication  Systems,  Artech 
(1981). 


The  H  matrix  Cor  N  =  3  is 


H  = 


[: 


(97) 


Let  4>oi  denote  the  bearing  angle  from  station  i  at  coordinates  (x^y^)  to  the 
reference  point  at  coordinates  (xq/Vq),  as  illustrated  in  figure  5.  Thus, 


*0i 


=  tan 


-1 


,  i  =  1,  2,  3 


(98) 
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Figure  5.  Angle  definitions  for  reference  and 
three  stations. 


Equation  ( 78)  may  be  expressed  as 


cos 

*01 

sin 

*01 
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cos 
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(99) 


The  covariance  matrix  P  can  be  evaluated  by  substituting  equations  (85), 
(96),  (97),  and  (99)  into  equation  (87).  The  components  of  P  defined  by 

equation  (53)  are 


23 


'A-u  '  '• 


(100) 


of  =  a[°t1  (sin  (j>02  -  sin  <t>03)2  +  a|2(sin  <t>01  "  sin  ^03)2 
+  a£3(sin  $01  -  sin  $02)2]  » 


°2  =  a[°tl(cos  (J>02  -  cos  <J>03)2  +  a2 2 (cos  $Q1  -  cos  $03)2 
+  o23(cos  $Q1  -  cos  4>02)2]  r 


(  101  ) 


°12  =  “[°tJcos  ^03  “  cos  ^02)(sin  *02  "  Sin  *03^ 

+  o22(cos  <J>03  -  cos  $oi  )(sin  $01  “  sin  *03) 

+  o23(cos  $02  -  cos  $01)(sin  $Q1  “  sin  *02)]  ' 


(102) 


where 


=  c2 [ ( cos  $01  -  cos  $02) (sin  $02  ■  sin  ^03^ 
-  (cos  $02  -  cos  $03) (sin  $Q1  -  sin  $02)] 


(103) 


If  any  two  bearing  angles  are  equal,  then  a2,  a|»  anti  °i  2  +  °°*  These  events 
correspond  to  reference  points  that  lie  along  a  line  passing  through  two  of 
the  stations. 

The  least-squares  or  maximum-likelihood  estimator,  determined  from  equa¬ 
tion  (86),  is 


x  =  xQ  +  /^  [(t,  -  D01/c)(sin  $02  "  sin  $Q3) 
+  (t2  -  DQ2/c)(sin  $03  -  sin  $Q1) 

+  (fc3  “  D03/cHsin  ♦oi  ”  sin  *02)]  ' 

y  =  y0  +  /^  [(t,  -  D01/c)(cos  $03  -  cos  $02) 
+  (t2  -  D02/c)(cos  $Q1  -  cos  $03) 

+  (fc3  *  D03/cKcos  ^02  "  cos  *01))  * 


(104) 


(105) 


To  determine  the  transmitter  range,  which  may  be  defined  as  the  distance 
between  the  transmitter  and  the  master  station,  it  is  convenient  to  align  the 
x-axis  with  the  line  j^etween  the  master  station  and  the  reference  point  and  to 
place  the  origin  of  the  coordinate  system  at  the  master  station.  If  the 
reference  point  is  near  the  transmitter  position,  then  x  is  a  suitable  range 


estimator  and  approximates  the  variance  of  the  range  estimator;  otherwise, 
the  range  can  be  estimated  by  (*2  +  y2)l/2,  a  suitable  estimator  for  the 


bearing  with  respect  to  the  x-axis  is 


=  tan 


The  estimator  bias  can  be  determined  from  equation  (15). 
linearization  error  and  using  equation  (83),  we  obtain 


(106) 


Neglecting  the 


bl  =  ^  (Efei](sin  4*o 2  ~  sin  *03)  +  E[e2](sin  *03 


sin 


*01  ) 


(107) 


+  E [ c 3 ] (sin  4>01  -  sin  4>0-> ) } 


b2  =  /a  {e[e1](cos  <j>03  -  cos  4,Q2)  +  e[e2](cos  <J>01  -  cos  <t>03) 
+  E [ e 3 ] (cos  4>02  -  cos  4>01)l  . 


(108) 


Nonzero  values  of  the  E[ejJ  are  caused  primarily  by  uncertainties  in  the 
station  positions,  synchronization  errors,  and  the  temperature  dependence  of 
the  receiver  delays  and  filter  characteristics. 

Assuming  that  n  =  He  has  a  Gaussian  distribution,  equations  (54),  (55), 

(73),  and  (100)  to  (103)  give  the  CEP  in  terms  of  the  bearing  angles  and  the 
arrival-time  variances.  For  a  fixed  deployment  of  stations,  the  locus  of 
transmitter  positions  with  a  constant  value  of  the  CEP  can  be  determined 
numerically.  For  this  purpose,  the  equations  may  be  expressed  in  terms  of  the 
Cartesian  coordinates  by  using  equation  (98)  and  it  is  assumed  that  the  refer¬ 
ence  point  coincides  with  the  transmitter  position  with  negligible  error  so 
that  Doi  =  Di. 

Let  L  denote  the  length  of  a  linear  array  of  three  stations  with  coordi¬ 
nates  (0,-L/2),  (0,0),  and  (0,L/2).  Assuming  that  the  lower  bound  of  inequal¬ 
ity  (91)  is  nearly  achieved  and  using  equation  (95),  we  obtain 

'n  exp[a(D0.  -  L)]  ,  i  =  1,  2,  3  ,  (109) 
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where  ofcL  denotes  the  lower  bound  of  0^  when  D0i  =  L*  It  is  assumed  that  a, 
n,  KE,  and  hence  are  identical  for  all  three  stations.  We  assume  that  the 

transmitter  and  the  stations  have  omnidirectional  antennas  so  that  KE  does  not 
depend  upon  the  bearing  angle  to  the  transmitter.  Figures  6  and  7  depict  loci 
of  constant  values  of  CEP/c 0^  for  a  =  0.  Only  the  first  quadrant  is  dis¬ 
played  because  of  the  symmetry  of  the  loci.  Figure  6  assumes  n  =  2,  which 
corresponds  to  free-space  propagation.  Figure  7  assumes  n  =  4,  which  might 
model  vhf  propagation  near  the  earth's  surface. 

In  figure  8,  the  stations  form  a  nonlinear  array  with  coordinates 
(0,-L/2),  (-L/2,0),  and  (0,L/2),  respectively.  The  most  significant  features 
are  the  singularities  in  the  values  of  the  CEP  along  the  lines  passing  through 
two  of  the  stations.  Consequently,  only  a  slight  spatial  nonlinearity  is 
permissible  if  a  broad  field  of  view  is  required.  However,  other  important 
factors  in  the  choice  of  station  positions  are  the  needs  to  maintain  line-of- 
sight  paths  from  potential  transmitter  positions  and  to  minimize  the  potential 
multipath  interference. 


Figure  6.  Loci  of  constant  CEP/ca, 
array  of  three  stations  with  n  = 


for  linear 


CEP/co. 


7.  Loci  of  constant  CEP/co 
of  three  stations  with  n  = 


Loci  of  constant  CEP/cotL  for 
array  of  three  stations  with  n 


Figure  8. 
nonlinear 
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In  figures  6  to  8,  the  parts  of  the  loci  for  which  X2  <  0.01 X1  are  indi¬ 
cated  by  dotted  lines.  For  these  small  values  of  X2/X1  t*'e  CEP  a  <lues“ 
tionable  measure  of  performance  of  the  passive  location  system.  A  more  suit¬ 
able  measure  may  be  the  length  of  the  major  axis  of  the  concentration  ellipse. 


2/ <X^  . 


It  follows  from  equation  (73)  that 


0.563L 


x2  <  o.oix1 


(ill  ) 


(112) 


where  <  is  given  by  equation  (60).  Thus,  the  dotted  lines  approximate  the 
loci  of  constant  values  of  Le/3.552/xc 

5.  LOCATION  USING  BEARING  MEASUREMENTS 

The  bearing  measurements  of  passive  direction-finding  systems  at  two  or 
more  stations  can  be  combined  by  a  direction-finding  location  system  to  pro¬ 
duce  an  estimate  of  transmitter  position.  The  transmitted  signal  may  be 
received  at  a  station  by  line-of -sight  propagation  or  after  atmospheric  re¬ 
flection  at  a  known  altitude.  A  single  bearing  angle  may  be  measured  at  each 
station  of  the  location  system.  Alternatively,  separate  azimuth  and  elevation 
angle  measurements,  possibly  made  by  orthogonal  interferometers,  can  be  used 
to  determine  transmitter  position.  In  the  absence  of  noise  and  interference, 
bearing  lines  from  two  or  more  stations  will  intersect  to  determine  a  unique 
location.  In  the  presence  of  noise,  more  than  two  bearing  lines  will  not 
intersect  at  a  single  point,  as  illustrated  for  a  planar  configuration  in 
figure  9.  Consequently,  processing  is  required  to  determine  the  optimal 
position  estimate.  Let  0i  denote  the  bearing  angle  measured  at  station  i 
relative  to  a  base  line  in  a  three-dimensional  coordinate  system  defined  so 
that  the  x-axis  is  parallel  to  the  base  line,  as  shown  in  figure  10.  If  the 
coordinates  of  the  station  are  (xiryj/zi)  and  the  coordinates  of  the  transmit¬ 
ter  are  (xt,yt,zfc),  then  in  the  absence  of  measurement  errors,  line-of-sight 
propagation  implies  that 


0^  =  cos 


■i  r  xt  -  xi 

V(xt  -  xi)2  +  (yt  -  yj2  +  (zt  -  zi)2 


,  0  <  0,  <  it  . 


(113) 


In  figure  10,  the  azimuth  angle  is  defined  in  the  plane  passing  through  the 
transmitter  and  perpendicular  to  the  z-axis.  It  is  positive  in  the  counter¬ 
clockwise  direction  relative  to  the  positive  x-axis.  If  the  elevation  angle 
of  the  station  relative  to  the  transmitter  is  known  approximately  or  is 
estimated  by  a  suitable  means,  such  as  a  vertical  interferometer,  then  ^  may 
be  calculated  using  the  geometrical  relation 


i 


COS  0^  =«  COS  <f>^  COS  , 

which  is  easily  derived  from  figure  10.  If  is  sufficiently  small, 
measured  bearing  is  well  approximated  by  the  azimuth,  which  is  defined  by 
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Figure  9.  Bearing  lines  from  three  direction¬ 
finding  systems. 
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Figure  10.  Angle  definitions  for  direction 
finding  system. 
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In  most  applications,  the  transmitter  is  known  to  lie  on  the  surface  of  the 
earth  or  at  a  fixed  altitude  so  that  zfc  is  known  and  does  not  have  to  be 
estimated.  Equation  (115)  is  used  in  the  estimation  of  the  (xt,yt).  The  use 
of  this  equation  is  equivalent  to  the  representation  of  the  three-dimensional 
problem  by  a  two-dimensional  model.  In  the  model,  the  transmitter  and  the 
stations  are  assumed  to  lie  in  the  same  plane  so  that  the  azimuths  are  identi¬ 
cal  to  the  bearings.  If  the  transmitter  and  the  stations  actually  lie  on  the 
earth's  surface,  the  model  is  an  idealization  that  neglects  the  curvature  of 
the  surface.  Two-dimensional  position  estimation  using  bearing  information  is 
often  called  triangulation. 

We  consider  in  detail  the  estimation  c r  two-dimensional  column  vector 

R  having  components  x  and  y.  Line-of-slght  propagation  is  assumed.  The 
measured  bearing  angle  4^  and  the  measurement  error  n^  satisfy 


c  fj_  (R)  +  n^  ,  i  =  1 ,  2, 


where 


f^(R)  =  tan' 


-  («i)  ■ 


i  *  1 ,  2,  .  .  .,  N  , 


and  the  station  coordinates  are  x^  and  y^.  In  matrix  form,  we  have 

+  -  f  ( R)  +  n  . 


(116) 


(117) 


(118) 


Let  the  column  vector  Rg  with  components  xQ  and  yQ  specify  a  reference  point, 
which  may  be  chosen  to  be  in  the  middle  of  the  polygon  bounded  by  the  measured 
bearing  lines.  Let  4>Qi  denote  the  bearing  angle  from  station  i  to  the  refer¬ 
ence  point.  Then, 


•  a  yp  -  yj 

sin  4>oi  =  Doi 


x0  ~  xi 


9  i  s  1»  2f  *  »  »>  N  t 


where 


D0i  -  [(xo  -  xi)2  +  (yo  -  yi)2]l/2  »  i  -  i,  2, 

equation  (7)  with  x  =  R  and  Xg  =  Rg,  we  obtain 


-(sin  4>01)/D01  (cos  ♦01)/D01 


.,  N  . 


(119) 


(120) 


(121  ) 


(sin  *on)/Don  (cos  ♦0h)/Dc 


The  least-squares  or  maximum-likelihood  estimator  is 

*  a  Rp  +  (GTH"lG)-l  GTR-l*r  - 


(122) 
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where  H  is  the  covariance  matrix  of  the  bearing  measurement  errors  and 


♦r  =  ♦  "  £(*o)  * 


The  ith  component  of  4r  is 


i  =  1,  2,  .  .  .,  N 


(123) 


(124) 


which  is  the  bearing  angle  relative  to  the  line  between  station  i  and  the 
reference  point,  as  depicted  in  figure  11. 
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Figure  11.  Geometry  of  transmitter, 
reference  point,  and  a  station. 


If  the  bearing  measurement  errors  are  independent  random  variables  with 
variances  o^,  i  =  1 ,  2,  .  .  .,  N,  then 


(125) 


Direct  calculation  using  equations  (16),  (121),  and  (125)  establishes  that  the 
elements  of  the  covariance  matrix  of  R  are 

o?  *  E[(x  -  x)2]  *  - - — ~  ,  (126) 

1  1  1  v\  -  v2 

ai  »  E[(y  -  y)2] - -  ,  (127) 

yX  - 

o,2  -  E[(x  -  x)(y  -  y)]  -  —  y_  —  ,  (128) 


where 
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It  follows  from  equations  (121),  (122),  and  (125)  that  the  components  of 
the  linearized  least-squares  estimator  are 


1  y  (v  cos  ♦oi  -  V  sin  *01) 

"  x0  +  —  “7  A,  *ri  - 1  I? 


pX  -  v2  i=1  ri  Doic|i 


y  “  Yo 


1  ?  .  (X  cos  *Qi  -  v  sin  <t>oi) 

+  "T  T5  .1.  ♦ri  I - "5  * 


pX  -  v2  i=1  ri 


D0i°|i 


(132) 


(133) 


Similarly,  if  the  linearization  error  is  negligible,  the  bias  components  are 


J-j  .1  *Kl 

-  v*  1=1 


(v  cos  4>oi  -  p  sin  4>oi) 
D0i°|i 


1  !  ,fnl  (X  cos  ♦pi  -  v  sin  »Qi) 

“  iil  [niJ  Dn,o?, 
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(134) 
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The  dependence  of  the  estimator  and  bias  on  o^,  i  =  1,  2,  .  .  .,  N,  is 
eliminated  because  of  cancellation  in  equations  (132)  to  (135)  if  these  vari¬ 
ances  are  all  equal.  This  equality  is  a  reasonable  assumption  if  the  receiv¬ 
ers  are  identical  and  much  closer  to  each  other  than  to  the  transmitter. 

Let  p^  denote  the  shortest  distance  from  the  reference  point  to  the  meas¬ 
ured  bearing  line  at  station  i,  as  depicted  in  figure  11.  Suppose  that  the 
reference  point  is  close  to  the  true  transmitter  position,  and  that  the  meas¬ 
urement  errors  are  small.  Then 


1,  2,  .  .  . ,  N  , 


(136) 


cos  ■  cos  4^,  sin  ■  sin  ,  i  =  1,  2,  .  .  .,  N 


(137) 
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Substituting  equations  (136)  and  (137)  into  equations  (129)  to  (133),  we 
obtain  the  components  of  an  estimator  that  depends  upon  the  measurements  p^ 
and  i  =  1,  2,  .  .  N.  This  estimator,  called  the  Stansfield  algorithm, 

was  originally  derived  from  heuristic  arguments  and  the  assumption  of  small 
bearing  measurement  errors . ^  if  R^  is  close  to  R,  then  the  linearized  least- 
squares  estimator  is  preferable  to  the  Stansfield  algorithm,  which  produces  a 
larger  estimator  bias  unless  the  bearing  errors  are  small.  However,  if  the 
bearing  errors  are  large,  it  may  not  be  possible  to  choose  *0  close  to  R.  In 
this  case,  it  is  not  clear  which  estimator  is  preferable. 

The  mean-square  ranging  error  associated  with  direction-finding  systems  is 
defined  as  the  average  variance  of  Doi^ri: 


=  H  J,  D0i°|i 


(138) 


In  analogy  to  equation  (89),  the  GDOP  associated  with  an  unbiased  estimator 
and  a  direction-finding  location  system  is  defined  as 


(139) 


If  the  geometry  is  such  that  the  bearing  variances  are  nearly  equal,  then  the 
GDOP  is  only  weakly  dependent  on  them.  From  equation  (74),  it  follows  that 


CEP  »  (o.750j)gdoP  . 


(140) 


The  variance  of  a  bearing  estimator,  a?,  is  due  primarily  to  the  thermal 
and  environmental  noise.  Approximate  expressions  for  o?  are  known  for  various 
direction-finding  systems  operating  in  white  Gaussian  noise.6 7  In  most  cases, 
if  E/Nq  is  sufficiently  large,  o|  can  be  expressed  in  the  form 


-k*)-1  ■ 


(141 ) 


where  8?  is  a  function  of  the  system  parameters  other  than  E/Nq,  and  the 
variation  of  the  signal  energy  with  the  distance  to  the  transmitter  can  be 
determined  from  equation  (95).  For  example,  consider  a  planar  configuration 
and  a  phase  interferometer  with  its  antennas  pointing  in  the  direction  of  the 
positive  x-axis.  It  can  be  shown  that  if  the  estimator  bias  is  small,  then6 
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60.  J.  Torrieri,  Principles  of  Military  Communication  Systems,  Aztech 
(1981). 

7C.  J .  Ancker,  Airborne  Direction  Finding — The  Theory  of  Navigation  Errors, 
IRE  Trans .  Aeronaut.  Navig.  Electron.  ANE-5  (December  1958),  199. 


where  fg  is  the  carrier  frequency  of  the  received  signal,  d  is  the  maximum 
separation  between  the  interferometer  antennas,  and  $  is  the  true  bearing 
angle. 

As  a  specific  example,  we  consider  identical  stations  that  are  symmetri¬ 
cally  located  with  respect  to  the  reference  point  so  that 


=  i+1 )  '  i  **  1  #  2,  .  .  .,  int(N/2)  ,  (143) 

DOi°|i  “  °0(N-i+1 ) °|(N-i+1 )  •  1  2»  •  •  •»  int(N/2)  ,  (144) 

where  int(x)  denotes  the  largest  integer  in  x.  If  N  is  odd,  we  further  assume 
that 

(J>oi  =  0  ,  i  -  int(N/2)  +  1  ,  N  is  odd.  (145) 

A  possible  configuration  for  N  -  5  is  illustrated  in  figure  12.  This  example 
is  probably  unrealistic  for  ground  stations  if  N  >  4,  but  might  adequately 
represent  location  estimation  by  an  aircraft  that  samples  bearing  data  at 
evenly  spaced  points  along  its  trajectory.  Substitution  of  equations  (143), 
(144),  and  (145)  into  equation  (131)  yields  v  =  0,  which  implies  that  o12  ” 
0.  We  conclude  that  the  symmetrical,  but  not  necessarily  linear,  placement  of 
the  stations  with  respect  to  an  accurately  located  reference  leads  to  uncor¬ 
related  coordinate  estimates.  For  an  aircraft,  o12  *  0  allows  interpretation 
of  o|  as  the  variance  of  the  "cross-range"  estimation  error  and  as  the 
variance  of  the  "down-range"  estimation  error. 


Figure  12.  Configuration  of  five  symmet 
rically  located  stations. 


If  N  «*  2,  equations  (126)  to  (133)  give 
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y  =  Vo  +  —cos  ♦Q-  (*t1  +  *r2) 


(149) 


If  the  reference  point  is  located  at  the  intersection  of  the  two  measured 
bearing  lines,  then  =  $r2  m  0.  It  follows  that  (x,y)  *  (xo>yo)'  as  ex“ 

pected.  From  equations  (138),  (139),  (146),  and  (147),  we  obtain 


GDOP 


n 

sin  2$q1 


(150) 


The  minimum  value  of  the  GDOP,  equal  to  /7,  is  attained  when  $qi  “  x/4.  Since 
01 2  “  0,  equations  (54),  (55),  and  (73)  give 

CEP  =  0.563  max(o1,o2)  +  0.614  minfo^Oj)  •  (151) 
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If  N  =  3,  the  variance  of  x  remains  the  same,  but  the  variance  of  y  becomes 
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which  shows  that  the  extra  station  only  improves  the  estimation  of  the  y- 
coordinate  of  the  transmitter.  As  the  transmitter  range  increases,  $Q1  de¬ 
creases  and  thus  o2/o^  increases . 

If  n  in  equation  (118)  has  a  Gaussian  distribution,  equations  (54),  (55), 
(73),  and  (126)  to  (131)  give  the  CEP  in  terms  of  the  bearing  angles  and  their 
variances.  Assuming  that  the  reference  point  coincides  with  the  transmitter 
position  so  that  -  D^  and  is  equal  to  the  bearing  angle  to  the  trans¬ 

mitter  position,  the  locus  of  positions  with  a  constant  CEP  can  be  determined 
numerically  by  using  equations  (119)  and  (120). 


Consider  a  linear  array  of  three  stations  with  coordinates  (0,-L/2), 
(0,0),  and  ( 0 , Ii/ 2 ) •  Each  station  has  an  interferometer  with  omnidirectional 
antennas  pointing  in  the  direction  of  the  positive  x-axis.  Let  0..  denote  the 
value  of  0^  when  -  L  and  ■  0.  Assuming  that  o.^  is  identical  for  all 

three  stations  and  that  the  lower  bound  of  inequality  (142)  is  nearly 
achieved,  equation  (95)  yields 
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O  STATION 

—  X2/X1<0.01 

—  Xtf  Xi  >0.01 


Figure  14.  Loci  of  constant  CEP/L0^L  for  linear 
array  of  three  stations  with  n  =  4. 


In  figure  15,  the  stations  form  a  nonlinear  array  with  coordinates 
(0,-L/2) ,  (-L/2 ,0 ) ,  and  (0,L/2).  A  comparison  with  figure  8  indicates  that 
the  adverse  effect  of  the  nonlinear  configuration  is  usually  less  for  direc¬ 
tion-finding  systems  than  for  hyperbolic  systems. 

Figure  16  plots  the  constant  CEP/L0aL  loci  for  a  linear  array  of  five 
stations  with  coordinates  (0,-L/2),  (0,-L/4),  (0,0),  (0,L/4),  and  (0,L/2).  A 
comparison  with  figure  14  shows  the  CEP  improvement  from  adding  two  stations 
while  maintaining  a  constant  baseline  length  equal  to  L.  In  general,  the  CEP 
is  roughly  inversely  proportional  to  /!T. 

For  two-dimensional  transmitter  location  with  three  stations,  a  comparison 
of  figures  13  to  15  with  figures  6  to  8  indicates  that  for  hyperbolic  systems 
to  offer  a  significant  performance  advantage  over  direction-finding  location 
systems  when  a  >*  0,  it  is  necessary  that 


qc0tL  <  L0^l 


(155) 


where  q  »  5.  Substituting  equations  (154)  and  (110)  and  assuming  equal  param¬ 
eter  values  for  the  two  systems,  we  obtain  the  criterion 


n  irqf-d  <  LB. 


(156) 


Figure  15.  Loci  of  constant  CEP/La..  for 
nonlinear  array  of  three  stations  with  n 


CEPfLo 


Figure  16.  Loci  of  constant  CEP/Lo^.  for  linear 
array  of  five  stations  with  n  =  4. 
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Consequently,  for  the  radar  signal  leading  to  equation  (93),  hyperbolic 
systems  offer  a  potential  advantage  only  if 


Tp(*qf0d)2  <  BL2  .  (157) 

For  the  communications  signal  leading  to  equation  (94),  a  significant  advan¬ 
tage  requires 


SZ  qfQd  <  BL 


(158) 


Inequalities  (157)  and  (158)  indicate  that  hyperbolic  systems  increase  in 
desirability  as  the  array  length  and  signal  bandwidth  increase. 

6.  OTHER  LOCATION  METHODS 


When  the  receivers  are  moving,  it  may  be  possible  to  use  the  known  receiv¬ 
er  trajectories  to  enhance  the  accuracy  of  the  transmitter  location.  For 
example,  three  bearing  measurements  and  two  turns  by  an  aircraft  can  be  used 
to  greatly  reduce  the  effect  of  strong  unknown  biases  in  the  measurements.® 

Moving  receivers  can  exploit  the  Doppler  shift  in  several  ways.  In  the 
absence  of  noise,  the  measured  frequency  at  a  receiver,  fm,  is  related  to  the 
transmitted  frequency,  ffc,  by 


f 


m 


+ 


ffcv  cos  $ 
c 


(159) 


where  c  is  the  signal  velocity,  vf  is  the  velocity  of  the  receiver  in  the 
direction  to  the  transmitter,  v  is  the  receiver  velocity,  and  <J>  is  the  bearing 
angle  to  the  transmitter  relative  to  the  velocity  vector,  as  shown  in  figure 
17.  Therefore,  the  bearing  angle  can  be  estimated  if  fm  is  measured  and  ffc, 
v,  and  c  are  known.  Bearing  measurements  from  several  receivers  can  be  com¬ 
bined  to  obtain  a  transmitter  location  estimate  as  in  section  5.  Another 
approach,  which  may  be  less  sensitive  to  inaccuracies  in  the  assumed  value  of 
ft,  is  to  measure  the  Doppler  difference,  which  is  defined  as 

fm1  “  Cm2  *  c^  (V1  cos  *1  -  v2  cos  *2)  '  (160) 

where  the  subscripts  1  and  2  refer  to  receivers  1  and  2.  The  differential 

Doppler  is  defined  as  the  integral  of  fm1  -  fm2  over  time.  If  ft  does  not 

change  too  rapidly  over  the  integration  interval,  the  differential  Doppler  is 


f^2  (fm1  "  fm2)  dt  *  [°1  C^2)  "  ~  ViM  +  '  (161) 


®M.  Mangel,  Three  Bearing  Method  for  Passive  Triangulation  in  Systems  with 
Unknown  Deterministic  Biases,  IEEE  Trans.  Aerosp.  Electron.  Syst.  AES- 17 
(November  1981),  814. 


where  ffca  is  the  average  transmitted  frequency  and  D^(tj),  i,  j  =  1,  2,  is  the 
distance  of  receiver  i  from  the  transmitter  at  time  j.  The  right-hand  sides 
of  equations  (160)  and  (161)  can  be  expressed  in  terms  of  the  transmitter 
coordinates.  Thus,  in  the  absence  of  noise,  a  Doppler  difference  or  a  differ¬ 
ential  Doppler  measurement  determines  a  surface  on  which  the  transmitter  must 
lie.  A  location  estimator  can  be  derived  in  a  manner  analogous  to  the  deriva¬ 
tions  of  sections  4  and  5.  Because  of  the  need  for  a  precise  estimate' of  ffc 
or  f^#  Doppler  location  systems  appear  to  be  most  useful  in  the  location  of 
transmitters  of  narrowband  signals. 


TRANSMITTER 


Figure  17.  Moving  receiver. 


Doppler,  arrival-time,  and  bearing  measurements  at  the  same  or  different 
receivers  can  be  combined  in  hybrid  location  systems.  The  combined  measure¬ 
ments  may  allow  a  reduction  in  the  number  of  receivers  required  for  a  given 
location  accuracy  and  may  facilitate  the  resolution  of  ambiguities. 

To  accommodate  a  moving  transmitter,  the  observation  interval  can  be 
decreased  so  that  the  transmitter  is  nearly  stationary  during  the  interval  and 
points  on  the  trajectory  can  be  located.  However,  decreases  in  the  observa¬ 
tion  interval  eventually  lead  to  unacceptably  large  estimation  errors,  and 
other  methods  must  be  adopted.  If  the  trajectory  can  be  described  by  a  low- 
order  polynomial  in  time  and  if  a  sufficient  number  of  stations  or  measure¬ 
ments  are  available,  it  is  possible  to  estimate  the  coefficients  by  expanding 
the  dimension  of  the  estimator  x.  Alternatively,  if  the  differential  equa¬ 
tions  of  motion  are  known,  Kalman  filters  can  be  used  to  track  the  transmitter 
movement.9' 10  However,  the  implementation  complexity  of  a  passive  location 
system  with  Kalman  filters  is  usually  considerably  greater  than  that  of  a 
hyperbolic  or  direction-finding  location  system  for  stationary  transmitters. 


9A.  Gelb,  ed.,  Applied  Optimal  Estimation,  MIT  Press  (1974). 

10P.  S.  Maybeck,  Stochastic  Models,  Estimation,  art!  Control,  volume  1, 
Academic  Press  (1979). 
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